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ABSTRACT 

We investigate numerical solution of Q 2 evolution equations for structure functions 
in the nucleon and in nuclei. (Dokshitzer-Gribov-Lipatov-)Altarelli-Parisi and Mueller- 
Qiu evolution equations are solved in a brute-force method. Spin-independent flavor- 
nonsinglet and singlet equations with next-to-leading-order a s corrections are studied. 
Dividing the variables x and Q 2 into small steps, we simply solve the integrodifferential 
equations. Numerical results indicate that accuracy is better than 2% in the region 
10~ 4 < x < 0.8 if more than two-hundred Q 2 steps and more than one-thousand x 
steps are taken. The numerical solution is discussed in detail, and evolution results 
are compared with Q 2 dependent data in CDHSW, SLAC, BCDMS, EMC, NMC, 
Fermilab-E665, ZEUS, and HI experiments. We provide a FORTRAN program for 
Q 2 evolution (and "devolution") of nonsinglet-quark, singlet-quark, + g i; and gluon 
distributions (and corresponding structure functions) in the nucleon and in nuclei. This 
is a very useful program for studying spin-independent structure functions. 
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Program Summary 



Title of program: BF1 

Computer: AlphaServer 2100 4/200; Installation: The Research Center for Nuclear 
Physics in Osaka 

Operating system: Open VMS V6.1 
Programming language used: FORTRAN 77 
Peripherals used: Laser printer 

No. of lines in distributed program, including test data, etc.: 2439 

Keywords: Structure function, parton distribution, Q 2 evolution, numerical solution. 

Nature of physical problem 

This program solves Altarelli-Parisi equations or modified evolution equations (Mueller- 
Qiu) with or without next-to-leading-order a s effects for a spin-independent structure 
function or quark distribution. Both flavor- nonsinglet and singlet cases are provided, 
so that the distributions, xq NS , xq s , xqf = xqi + xqt (?=quark flavor), xg, xF NS , xF s , 
and xFf in the nucleon and in nuclei can be evolved. 

Method of solution 

We divide the variable x (and Q 2 ) into very small steps, and integration and differ- 
entiation are defined by = \-f( Xm + 1 } f( x ™)] an( j f d x f( x ) — ^ Ax m f(x m ). 

UX L\X m J m=l 

Then, the integrodifferential equations are simply solved step by step, and this method 
is so called brute-force method. If the step numbers are increased, accurate results 
should be obtained. 

Restrictions of the program 

This program is used for calculating Q 2 evolution of a spin-independent flavor-nonsinglet- 
quark, singlet-quark, qf , and gluon distributions (and corresponding structure func- 
tions) in the leading order or in the next-to-leading-order of a s . Q 2 evolution equations 
are the Altarelli-Parisi equations and the modified ones (Mueller-Qiu). The double pre- 
cision arithmetic is used. The renormalization scheme is the modified minimal subtrac- 
tion scheme (MS). A user provides the initial structure function or quark distribution 
as a subroutine or as a data file. Examples are explained in sections 4.2 and 4.3. Then, 
the user inputs twenty-one parameters in section 4.1. 

Typical running time 

Approximately five minutes on AlphaServer 2100 4/200 in the nonsinglet case, sixty 
minutes in the singlet-quark evolution, and eighty minutes in the singlet evolution with 
recombination effects. 
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LONG WRITE-UP 



1. Introduction 

Subnucleon degrees of freedom could be investigated in high-energy lepton-nucleon 
interactions. Measured structure functions depend in general on two variables, Q 2 = 
—q 2 and x = Q 2 /2P-q, where q is the four-momentum transfer and P is the nucleon mo- 
mentum. Bjorken scaling hypothesis suggests that structure functions are independent 
of Q 2 . However, the scaling is not exactly satisfied and they have weak logarithmic 
Q 2 dependence, which is refer to as scaling violation. Although the structure func- 
tions themselves cannot be calculated except for lattice QCD methods, it is possible 
to estimate their Q 2 variations within perturbative QCD. 

An intuitive way of describing the scaling violation is to use the (Dokshitzer-Gribov- 
Lipatov-)Altarelli-Parisi equation [3J. For example, the flavor-nonsinglet Altarelli- 
Parisi (or DGLAP) equation is given by 



where q NS (x,Q 2 ) is a nonsinglet quark distribution, P NS (x) is a nonsinglet splitting 
function, and a s (Q 2 ) is the running coupling constant. The integrodifferential equation 
describes the progress that a quark with the nucleon's momentum fraction y radiates 
a gluon and it becomes a quark with the momentum fraction x. The splitting function 
determines the probability of the splitting process. The singlet Altarelli-Parisi equa- 
tions are more complicated due to gluon participation in the evolution process. The 
flavor-singlet part is given by coupled integrodifferential equations, which are discussed 
in section 2. 

We study not only the Altarelli-Parisi equations but also modified evolution equa- 
tions due to parton recombinations (Mueller and Qiu) 0. The recombination mech- 
anism produces additional terms in the Altarelli-Parisi evolution equations. Because 
gluons play a major role at small x, gluon recombination processes are investigated in 
Ref. ||. Although the nonsinglet equation remains unchanged, the singlet equations 
are modified due to the gluon recombinations. The additional terms are nonlinear, so 
that it is not obvious how to solve the evolution equation numerically. On the other 
hand, the parton recombination mechanism is increasingly important due to recent 
HERA data in the small x region. It is interesting to test whether the F 2 data at small 
x could be related to higher-twist effects such as the recombination contributions 0. 
In addition, the recombination mechanism is used for explaining nuclear shadowing |4| , 
and it is an important factor in Q 2 evolution of nuclear structure functions. 

The Q 2 evolution equations are often used in experimental analysis and also in the- 
oretical calculations, so that it is worth while creating a computer program in solving 
the equations accurately. A number of methods have been developed and they 




(1.1) 
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include brute-force methods, Mellin-transformation methods, and orthogonal polyno- 
mial methods. Among these methods, the Laguerre-polynomial one || is considered an 
efficient one in computing time and in numerical accuracy. We have been investigating 
the Laguerre method and find that accuracy may not be very satisfactory at small x 
and at large x particularly in the nonsinglet case. With development of high-energy 
accelerators, the small x region becomes important. In fact, the x region of 10 -4 could 
be reached with large enough Q 2 at which perturbative QCD could be used. Consid- 
ering structure-function studies at HERA and at future high-energy laboratories, we 
should have a computer program which enables the evolution with good accuracy at 
small x. In addition to this problem, there is other difficulty in the Laguerre method. 
It is not obvious how to apply the Laguerre method to the case of modified Q 2 evolu- 
tion equations with parton recombinations. It is because of the existence of nonlinear 
recombination terms in the evolution equations. 

These issues motivated us to explore an alternative method. As a possible way 
to solve the above difficulties, we study a brute-force method. This method was, for 
example, investigated in Ref. f7|; however, there is no published article in discussing 
the details of the solution and its accuracy. The method is perhaps the simplest 
one in solving the integrodifferential equations. We divide the variables x and Q 2 

into very small steps, and integration and differentiation are defined by df(x)/dx = 

, n x 

[f(x m+ i) — f(x m )]/Ax m and / dxf(x) = ^ Ax m f(x m ). Then, the evolution equation 

m=l 

can be solved step by step. If N x is large enough, we should be able to get accurate 
numerical results. Our research purposes are to investigate the details of numerical 
accuracy and to provide a useful computer program in the brute-force method. 

The evolution equations, which are solved by the brute-force method, are given in 
section 2. We explain the brute-force method in section 3. In section 4, information 
for running the program BF1 is supplied. Each subroutine in the program is explained 
in section 5. Numerical results and comparisons with experimental data are given in 
section 6. The results are summarized in section 7. Explicit equations of the splitting 
functions and other necessary quantities are given in appendices. 
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2. Q 2 evolution equations 



We solve two types of evolution equations. One is the ordinary equations, so called 
Altarelli and Parisi equations, and the other is the modified ones due to recombinations 
by Mueller and Qiu. Parton distributions or structure functions in the nucleon can be 
evolved in both methods. The modified equations can also handle evolution of nuclear 
parton distributions. 



2.1 (Dokshitzer-Gribov-Lipatov-)Altarelli-Parisi equations 



The nonsinglet "Altarelli-Parisi" equation is given in Eq. (1.1), which is valid both 
in the leading order (LO) case and in the next-to-leading order (NLO) one. NLO effects 
can be included in the running coupling constant a s (t) and in the splitting functions 
Pij(z). In order to remove the extra Q 2 dependence in front of the integral in Eq. (1.1), 
we use the variable t defined by 



a s (Q 2 ) 
a s (Q 2 ) 



(2.1) 



instead of Q 2 . The parton distribution and the splitting function multiplied by x satisfy 
the same integrodifferential equation. Therefore, defining f(x) by 



f{x) = xf(x) 



we rewrite the evolution equation as 





Of 



q NS [x, t) 



dy 

x y 



(2.2) 



(2.3) 



The singlet case is more complicated than the nonsinglet one due to the gluon 
participation in the evolution, and it becomes coupled integrodifferential equations. 
The singlet quark distribution is defined by q s (x, t) = q^(x, t), where i is the quark 

i 

flavor and the qf distribution is given by q^(x, t) = qi(x, t)+q i {x 1 1). We write evolution 
equations in terms of the qf distribution R H|: 



|^ + (x,t) 



d_ 

Of 



g{x,t) 



l dy 

x y 

x y 



<J<I, 



X 



g, + (y,t) + P. 



x 



99 



g(y,t) 



[2Aa) 



(2.46) 



Each term in Eqs. (2.4a) and (2.4b) describes the process that a parton pj with the 
nucleon's momentum fraction y splits into a parton pi with the momentum fraction 
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x and another parton. The splitting function P PiP] (z) determines the probability that 
such a splitting process occurs and the p^-parton momentum is reduced by the fraction 



z. Using P 
we obtain 
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9 fat) 
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p, 



9<r 



Q s (y,t) + Pc 



X 



99 



g(y,t) 



(2.5a) 
(2.56) 



x 

All the splitting functions in Eqs. (2.5a) and (2.5b) are listed in Appendix B. If the 
summation over % is taken in Eq. (2.5a), it becomes the evolution equation for the 
singlet distribution q s : 



9 ~ , ^ 

ol % fa t) 



x y 



Pn 



NS 



- + 2N f C F T R F qq ( - 




+ ivyp ( 



X 



f r d9 



- g(v,t) 
yj 



(2.6) 



The reason for writing the evolution equations in term of the flavor- dependent distri- 
bution qf instead of the singlet one q s is because antiquark distributions are neither 
SU(3) flavor nor SU(2) flavor symmetric. 

Next, we discuss NLO effects in the evolution equations. NLO effects are included 
in the coupling constant oc s (t), in the splitting functions P p . p .(z), and in coefficient 
functions. Explicit NLO expressions are given in appendices. It should be noted 
that the LO or NLO expression of a s (Q 2 ) is used in the LO or NLO evolution case 
respectively in Eq. (2.1). Once the NLO corrections are included in the evolution, the 
renormalization scheme has to be specified. Throughout this paper, we use the MS 
scheme. The splitting function is given by the LO and NLO ones as 



Then the splitting functions Pij(z) in Eqs. (2.4a) and (2.4b) are 

P i ,( x )=pM( x ) + ^& Rij ( x ) , 



2tt 



where the function Rij(z) is 



ft 5(0), 



R ij {x)^P^\x)-^P^{x) . 



(2.7) 



(2.8) 



(2.9) 
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The second term in Eq. (2.9) appears because of the transformation from Q 2 to t. To 
be precise, the splitting functions P^ in Eqs. (2.3)-(2.6) and (2.8) should be denoted, 



However, we 



for example, P[a because it is different from xPij = xP^ +(a s /2jr)s^ 
omit the prime throughout this paper for using simpler notations. 

We discussed the evolution equations for quark and gluon distributions. In calcu- 
lating structure functions, a correction due to the NLO coefficient function should be 
taken into account. In the nonsinglet case, it is given by 



iNS, 



x,Q 2 



dy 

y 



X 



M 



a. 



NS (y,Q 2 ) 



[2.10) 



where n denotes the type of a structure function, F n =xFi, F2, or xF% for 71= 1, 2, or 
3 respectively, and C% is a quark coefficient function. In the singlet case or in the qf 
case, an additional gluon correction term should be taken into account: 



dy 



X 



,®s qt(y,Q 2 ) + 



dy 



x 



a. 



g(y,Q 2 



(2.11) 



where is a gluon coefficient function. and C® are given in Appendix C. The 
Q 2 evolution of a quark distribution is first calculated, then the structure function at 
Q 2 is evaluated by using the convolution integrals in Eq. (2.10) or (2.11), so that this 
method is called "two-step evolution". 

The above procedure requires an initial quark (and gluon) distribution for getting 
a structure function at certain Q 2 . It cannot be used for the NLO evolution if a struc- 
ture function is given as the initial distribution. In such a case, "one-step-evolution" 
equations for the structure function are useful. They are derived from the evolution 
equations and the convolution equations with the coefficient functions. The nonsinglet 
equation is in the same form with Eqs. (2.3) and (2.8) M: 



— r NS {x, t) 



x dy_ 

y 



5(0) ( £ 

NS\ y 



+ [y 



F N3 (y,t) 



except for taking 



(2.12) 



(2.13) 



B^(x) (ti=1, 2, or 3) is the a s correction in the coefficient function in Appendix C. 
Using the above one-step equation, we can evolve the nonsinglet structure function 
itself. In the similar way, the one-step F^~ evolution equation is given by H] 



with 



F, + (*,*) 
g(x,t) , 

R(x) 



'dy 

y 



p(0) 



+ 



a s (t) 



2tt 



R 



\ + (y,t) 
g{y,t) 



p(D, 



X 



fj D«(x) + E(x) 



(2.14) 



(2.15) 
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P = P/x, = D^/x, and E = E/x are given by 



p(0),(l) 



(0),(1) 

git 



(x) P$M{x 



J 



and 



D (D 



B*(x) Bi{x) 




#99 (^) 

Egg (X) 



(2.16) 



(2-17) 



(2.18) 



Each matrix element in Eqs. (2.16), (2.17), and (2.18) is given in Appendix B. It should 
be mentioned that the above one-step-evolution equation is slightly different from the 
one in Ref. ||. This discrepancy is because we use a convention for the momentum 

conservation, / dx\S^x{qi + qA + xg] = 1, which is different from the one used by 

Herrod, Wada, and Webber, / dx[Fs + xg'] = 1. Therefore, the matrices D < - 1 - ) and E 

Jo 

are different from those listed in their paper. 



2.2 Modified evolution equations due to recombinations (Mueller- Qiu equations) 

The Altarelli-Parisi evolution equations have been used extensively and they have 
been successful in describing many experimental data. However, as it becomes possible 
to reach the small x region by high-energy accelerators, it is necessary to investigate the 
details of small x physics. The gluon distribution g(x, Q 2 ) is a measurement of a gluon 
with transverse size ~ l/y/Q 2 , and the number of gluons per unit rapidity is xg(x, Q 2 ). 
If xg(x,Q 2 )/Q 2 « Rq (i?o=nucleon size) is satisfied, gluon interactions are neglected 
and the Altarelli-Parisi equations in Eqs. (2.3), (2.4), (2.5), and (2.6) can be used for 
the Q 2 evolution. However, if the transverse area xg(x, Q 2 )/Q 2 exceed the nucleon size 
~ Rq, gluons within a unit of rapidity spatially overlap, and the gluons are no longer 
considered as free partons. If the gluon density becomes large xg(x,Q 2 ) > RqQ 2 , the 
gluon interactions should be taken into account. These interactions are called parton 
recombinations, which are also used for explaining nuclear shadowing M. There are a 
number of studies on the recombinations. Among them, we employ the evolution equa- 
tions proposed by Mueller and Qiu [Q. They investigated gluon-gluon recombination 
effects on the evolution, and they proposed the following modified evolution equations: 
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i "':'- 41) °M 

) q s (y,t) + P, 



(2.196) 



gg r) g(y,t) 



T —a s (t)C 2 G f 1 ^ [g(y,t)] 2 \ , (2.19c) 
- 1 Jx y I 



*Ro 2 Q 2 I AT2 

where Ght(x, t) is a higher- dimensional distribution, which follows the evolution equa- 
tion: 

^G HT (x,t) = a s (t) Cl £j[9(y,t)] 2 ■ (2.19d) 



In order to solve this equation, the initial Gut distribution has to be supplied. Because 
Get is associated with the two gluon interactions, we may assume 



G HT (x,Ql) = K HT [g(x,Q 2 )] 2 



(2.20) 



where K HT is a constant. In Eq. (2.19a, b, c), R is the nucleon size, Nf is the 
number of flavor, N c is the number of color, and K—9/8 so that nuclear correction 

terms vanish for A = 1 [See Eqs. (2.22a, b, c)]. The splitting function P qg (z) is given 
in Appendix B. The recombination contributions enter into the evolution in the same 
a s order as the next-to-leading order in the splitting function, so that the NLO effects 
in the Altarelli-Parisi evolution terms should be included. 
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We comment on additional terms in Eq. (2.19). For the details the recombination 
processes, we refer the reader to the original paper B. The recombination terms are 
related to the two-gluon density, G^ 2 \x, Q 2 ) = G^'(xi,X2,Pit,P2t) with x\ = x 2 = x 
and p\ T = p\ T = Q 2 . Because the normalization of the two-gluon density is given 
by G ( - 2 \xi,X2,pl T ,P2 T )= (3R n /2)g(xi,pl T )g(x2,P2 T ) with n = it becomes 

G^(x,Q 2 ) = (3R n /2)[g(x, Q 2 )] 2 , which are the [g(x,Q 2 )] 2 terms in Eq. (2.19). It re- 
sults in the recombination factor, which is proportional to R^n^Ois/ Q 2 = 3a s / (AttRqQ 2 ) . 
The factor a s /Q 2 arises because a parton-parton fusion cross section is proportional to 
a s /Q 2 . This extra 1/Q 2 dependence is an interesting higher-twist effect, which could 
be investigated experimentally at small x and at small Q 2 . The 1/Rq factor arises due 
to the integration over the intrinsic transverse momentum. 

In the case of nuclear parton distributions, there are additional terms due to parton 
recombinations from different nucleons. Parton distributions in a nucleus are defined 
by 

p A (x, t) = p(x, t) - 5p A (x } t) , (2.21) 

where nuclear distributions are expressed by those per nucleon. The evolution of the 
first term is given in Eq. (19). The evolution of the nuclear correction terms is obtained 
by solving 
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(2.226) 
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Sg (x,t) 
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1 dy ~ fx 




) 



99 
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■kRq 



K 



2 



) 



Q 2 \n*-1 



1 / 87T 2 



<Xs{t) C\ 



(2.22c) 



where n is the number density n = 3A/(4irR\). Nuclear radius is given by Ra = R\A 1 ^ 
with Ri = 1.1 fm. For simplicity, we take R —Ri—1 fm and K—9/S so that the nuclear 
correction terms in Eqs. (2.22a, b, c) vanish at A = 1. We note that two parton number 
density is T^fai, X2, Q 2 ) = (3/2)RAnpi(xi)p2(x2) ■ Therefore, the overall A dependence 
is Rau oc A 1 / 3 , which is just the number of nucleons in the longitudinal direction. 
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3. Brute-force method 



It is useful to have a computer program of solving the evolution equations accurately 
because they are frequently used in theoretical and experimental studies. There are a 
number of methods such as Mellin-transformation and orthogonal-polynomial methods 
[Q. We have been studying a Laguerre-polynomial method and it is very efficient 
by considering computing time and numerical accuracy ||. However, the results are 
slightly worse in the nonsinglet case, particularly at small and large x. Furthermore, 
the evolution equations with the recombinations have complex nonlinear terms, which 
cannot be handled properly by the Laguerre polynomials. In addition to these issues, 
parton distributions with singular behavior at small x (xp(x) — > oo as x — > 0) could not 
be evolved either in the program LAG1 nor in LAG2NS ||. The singular behavior in 
sea-quark and gluon distributions attracts much attention recently due to the HERA 
F2 experimental data. 

Taking these difficulties into consideration, we decide to employ a brute-force 
method 0. The variables x and t are divided into small steps and integration and differ- 

entiation are defined by df(x)/dx = [f(x m+ i) — f(x m )]/Ax m and / dxf(x) = ^2 Ax m 

m=l 

f{x m ). Then, the evolution equations could be solved rather easily. For example, Eq. 
(2.3) is written in the following form: 

q NS (x k ,t j+1 ) = q NS {x k ,tj) + Atj E ^ ^ s (^) q NS (x m ,tj) . (3.1) 

m=k Xm Um/ 

If the initial distribution q NS (x m , t,=i) is given, the distribution at t\ = At is calculated 
in the above equation. Repeating this step N t — 1 times, we obtain the final distribution 
at t/Vf This method is very simple but N x and N t have to be large enough to get 
accurate results. However, we can reasonably expect that Nt does not have to be 
very large. This is because the scaling of structure functions works approximately, 
and they depend on the variable t weakly. On the other hand, N x has to be fairly 
large. Numerical problems are expected at small x if the x step is taken in the linear 
scale (Ax = 1/N X ). The small x region becomes increasingly important with the 
development of high-energy accelerators such as HERA. So it is necessary to have a 
good numerical method at small x as small as 10~ 4 . In order to satisfy this condition, 
the logarithmic-a; step A(log\ox) = \logiox m i n \/N x is taken in our analysis. 

The singlet and qf evolution equations in Eqs. (2.5a), (2.5b), and (2.6) are more 
complex, but they can be solved in the similar way. These equations are written in the 
brute-force method as: 
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Xk 



g{x m , tj) 

Qs \ X mi tj) 
g(yX m , tj) 



(3.26) 



(3.2c) 



If the initial distributions, q~i~(x m ,tj=o), q s (x m ,tj =0 ), and g(x m ,tj=o) are provided, 
evolved distributions at ti = At are calculated in the above equations. Repeating 
this step N t — 1 times, we obtain the final evolved distributions. The other evolution 
equations in section 2 can be solved in the same way. 

We did not write explicitly the integrals associated with 1/(1 — x) + terms in Eq. 
(3.2). However, special care has to be taken in calculating the 1/(1 — x) + terms in the 
splitting functions. The + function is defined in the integral region (0 < x < 1) [Eq. 
(B.2)], so that the integral is given in the brute-force method as 



Jx (I — X'). 



Y.Ax J^ J {1) + /(l)ln(l-s fc ) . 

m=k 



1 X m 

The integral with [ln(l — x)/(l — x)] + is given in the similar way: 



(3.3) 



dx'f(x') 



ln(l -x') 



x' 



N x 



= Ax m [f(x r , 

. + m=k 



1 Xm, £ 



(3-4) 

While the nonlinear recombination terms are difficult to be handled in the Laguerre 
method, it is straightforward to include them in the brute-force method. In addition, 
there is no evolution problem due to the singular behavior of parton distributions at 
small x, because all calculations in Eqs. (3.1) and (3.2) could be done in the x region, 
x > x m i n . x m i n is set for 10 -4 in our analysis, but it could be varied depending on one's 
interest. Because the recombination equations in section 2.2 are solved in the similar 
way, we do not write corresponding expressions in the brute-force method. 
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4. Description of input parameters 
and input distributions 



For running the FORTRAN-77 program BF1, a user should supply twenty-one input 
parameters from the file #10. In addition, an input distribution(s) should be given 
in a function subroutine(s) in the end of the FORTRAN program or in an input data 
file(s), #13, #14, #15, #16, and/or #17. Evolution results are written in the output 
file #11. We explain the input parameters and the input distributions in the following. 



4-1 Input parameters 



There are twenty-one input parameters. Numerical values of the parameters should 
be supplied in the file #10, then these are read in the main program. 



(1) IOUT 



(2) INPUT 

(3) IREAD 

(4) MODIFY 

(5) INDIST 

(6) IORDER 

(7) ITYPE 



x 



[or xF h 



x 



at fixed Q 2 (=Q2) in the file #11; 



(Q 2 ) [or xF NS (Q 2 )} at fixed x (=XX) in the file #11; 



1, write x and xq N: 

2, write Q 2 and xq 

3, x, xq s (x) [or xF s (x)], and xg(x); 

4, Q 2 , xq s (Q 2 ) [or xF s (Q 2 )], and xg(Q 2 ); 

5, x, xqf(x) [or xF^(x)}, xq s (x) [or xF s (x)], and xg(x); 

6, Q 2 , xqt{Q 2 ) [or xF+(Q% xq s (Q 2 ) [or xF s (Q 2 )], and xg(Q 2 ). 

1, parton distributions in the nucleon; 

2, parton distributions in a nucleus. 

Note: If MODIFY=2 and INPUT=2 are chosen, parton distributions 
both in a nucleus and in the nucleon should be supplied, 
give initial distribution(s) in function subroutine (s); 
read initial distribution(s) from data file(s). 
Altarelli-Parisi Q 2 evolution; 

Q 2 evolution with parton recombinations (Mueller-Qiu). 
do not write initial distribution(s); 
write initial distribution(s) in the file #12; 

write initial distribution(s) in the file #12 without calculating evolution, 
leading order in a s ; 
next-to-leading order, 
structure function xFx(x, Q 2 ); 

F 2 (x,Q 2 ); 
xF 3 (x,Q 2 ); 
quark distribution xq(x,Q 2 ). 
ITYPE determines the output distribution type. 



= 1 

= 2 

= 1 

= 2 

= 1 

= 2 

= 3 

= 1 

= 2 

= 1 

= 2 

= 3 

= 4 
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(8) ISTEP = 1, one-step evolution: xq{Ql) -> xq(Q 2 ) [or xF(Ql) -> xF(Q 2 )], 

^(Qo) 2#(Q 2 ); 

note: one-step xF^s evolution is not supplied; 
= 2, two-step evolution: xq(Q 2 ) ) — > xq(Q 2 ) — > xF(Q 2 ), 

xg{Ql) - x</(Q 2 ). 



(9) 


TA T / \ T > T ) 


= 1, q\ — q~i type distribution; 
= 2, q<i + type distribution. 


(10) 


Q02 


= initial Q 2 (= Qo in GeV 2 ) at which an initial distribution is supplied 






— id wnicn tne QistiiDUCion is evoiveQ ^ ^q)- 


(12) 


DLAM 


= QCD scale parameter Aq C d hi GeV. 


(13) 


NF 


= number of quark flavors (NF=3 or 4). 

= x at which Q 2 dependent distributions are written (IOUT=2, 4, or 6 


(14) 


XX 


(15) 


NX 


= number of x steps (NX<5000). 


(16) 


NT 


= number of t steps (NT<5000). 


(17) 


NSTEP 


= number of x steps or t steps for writing output distribution(s). 


(18) 


NXMIN 


= log w (minimum of x) [0 < min{x) = iq nxmin < XX\. 


(19) 


NA 


= mass number of a nucleus (NA=1 in the nucleon case). 


(20) 


RFM 


= #0=^1 (fm) in Eqs. (2.19) and (2.22). 


(21) 


DKHT 


= K HT in Ght(x) in Eq. (2.20). 
(Ght(x) is a higher dimensional gluon distribution). 



The meaning of IREAD is explained in section 4.2. IMORP=l or 2 means that the 
distribution is V] ca{qi + qi) or a>i{qi — qi) type, where are constants. For example, 

i i 

u v + d v —(u — u) + (d — d) and F% p + F^ p are q^ — q~i type distributions. F 2 ep is obtained 
by the convolution of the distributions, (A/9)x(u + u + c + c)+(l/9)x(d + d + s + s) and 
xg, with the corresponding coefficient functions, so that it is a qi + qi type distribution. 
In the evolution of nucleon structure functions, NA=1 should be supplied. In the 
Altarelli-Parisi evolution, the constants RKM and DKHT are unnecessary so that 
arbitrary constants may be supplied. 

For example, if one would like to evolve an initial singlet-quark distribution xq s at 
Q 2 =4 GeV 2 in the nucleon to the singlet structure function F 2) s at Q 2 =200 GeV 2 by 
the NLO Altarelli-Parisi with Nf—A and A=0.255 GeV, the input parameters could be 
IOUT=3, INPUT=1, IREAD=1, MODIFY=l, INDIST=1, IORDER=2, ITYPE=2, 
ISTEP=2, IMORP=2, Q02=4.0, Q2=200.0, DLAM=0.255, NF=4, XX=0.0, NX=1000, 
NT=200, NSTEP=100, NXMIN=-4, NA=1, RFM=0.0, and DKHT=0.0. In this case, 
the input file #10 is the following. 

3, 1,1,1,1, 2, 2, 2,2 

4.0, 200.0, 0.255, 4, 0.0, 1000, 200, 100, -4, 1, 0.0, 0.0 
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4-2 Input distributions supplied by function subroutines (IREAD=1) 



If IREAD=1 is chosen, an input distribution(s) at Qq should be supplied in the end 
of the FORTRAN program BF1 as a function subroutine (s). 

1) Nonsinglet case in the nucleon and nuclei 

An initial nonsinglet- quark distribution (or a nonsinglet structure function) at Q% 
should be given in QNSO(X) as a double precision function. As an example, the 
MRS(G) valence quark distribution xu v +xd v |L0| at Q§=4 GeV 2 is given in the original 
program BF1. 

2) Singlet case in the nucleon 

An initial singlet-quark distribution (or a singlet structure function) in the nucleon 
at Ql should be supplied in QSO(X), and an initial gluon distribution in the nucleon 
should be in G0(X). These subroutines are used in the nucleon case and also in the 
nuclear case with MODIFY=2. The MRS(G) +xd v +xS and xg distributions 

are given. 

3) qf (-F+) distribution case in the nucleon 

An initial qf distribution (or F+) in the nucleon at Ql should be supplied in QIO(X), 
a singlet-quark distribution (or a singlet structure function) in the nucleon should be 
in QSO(X), and an initial gluon distribution in the nucleon should be in G0(X). These 
subroutines are used in the nucleon case and also in the nuclear case with MODIFY=2. 
The MRS(G) distributions are given in each function subroutine, xqf = xd + xd is 
given as an example. 

4) Singlet case in a nucleus 

An initial nuclear singlet-quark distribution (or a singlet structure function) at 
Ql should be supplied in QSAO(X), and a nuclear gluon distribution should be in 
GAO(X). These subroutines are used only in the nuclear case. The SK xq^ a and xg Ca 
distributions Q are given as an example. 

5) qf (i^ + ) distribution case in a nucleus 

An initial nuclear qf distribution (or Ff) at Ql should be supplied in QIAO(X), 
a nuclear singlet-quark distribution (or a singlet structure function) should be in 
QSAO(X), and an nuclear gluon distribution should be in GAO(X). These subroutines 
are used only in the nuclear case. The SK distributions are given in each function 
subroutine. 
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4-3 Input distributions supplied by data files (IREAD=2) 



If IREAD=2 is chosen, an input distribution(s) at Ql should be supplied in a 
separate data file(s). 

1) Nonsinglet case in the nucleon and nuclei (data file #13) 

An initial nonsinglet-quark distribution (or a nonsinglet structure function) at Qq 
should be given in the data file #13 as shown in the following example. 



0.000100 
0.000110 
0.000120 
0.000132 
0.000145 



0.023868 
0.024942 
0.026069 
0.027251 
0.028491 



1.000000 



0.000000 



The first column is the x values and the second one is the corresponding xu v + xd v 
values. The data at x < x m i n and at £=1.0 must be supplied. 

2) Singlet case in the nucleon (data file #14) 

An initial singlet-quark distribution (or a singlet structure function) and a gluon 
distribution should be given in the data file #14 as shown in the following. 



0.000100 
0.000110 
0.000120 
0.000132 
0.000145 



3.137921 
3.114648 
3.091368 
3.068077 
3.044768 



23.163731 
22.485543 
21.825121 
21.181969 
20.555604 



1.000000 



0.000000 



0.000000 



The first column is the x values, the second is the xq s distribution, and the third 
is the xg distribution. The data at x < x min and at x=1.0 must be supplied. 

3) qf {Ff) distribution case in the nucleon (data file #15) 

An initial qf distribution (or Ff), a singlet-quark distribution (or a singlet structure 
function), and a gluon distribution should be given in the data file #15 as shown in 
the following. 
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0.000100 
0.000110 
0.000120 
0.000132 
0.000145 



1.229377 
1.220390 
1.211413 
1.202445 
1.193483 



1.000000 



0.000000 



3.137921 
3.114648 
3.091368 
3.068077 
3.044768 



23.163731 
22.485543 
21.825121 
21.181969 
20.555604 



0.000000 



0.000000 



The first column is the x values, the second is the xu + distribution, the third is 
xq s , and the fourth is xg. The data at x < x m i n and at x—1.0 must be supplied. 

4) Singlet case in a nucleus (data file #16) 

An initial nuclear singlet-quark distribution (or a singlet structure function) and a 
nuclear gluon distribution should be given in the data file #16 as shown in the nucleon 
case 2). In the MODIFY=2 case, the nucleon data file #14 should be also supplied. 

5) qf (Ff) distribution case in a nucleus (data file #17) 

An initial nuclear qf distribution (or Ff~), a nuclear singlet-quark distribution (or 
a singlet structure function), and a nuclear gluon distribution should be given in the 
data file #17 as shown in the nucleon case 3). In the MODIFY=2 case, the nucleon 
data file #15 should be also supplied. 
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5. Description of the program BF1 

5. 1 Main program BF1 

The main program reads twenty-one input parameters from the input file #10. The 
parameters are checked by the subroutine ERR. If there is an error, the program stops. 
Then, t defined in Eq. (2.1) is evaluated. Spline coefficients of the Spence function 
in Eqs. (B.7) and (B.8) are calculated by calling the SPLINE subroutine, and the 
interpolated function is used throughout the program for saving computation time. In 
the end, GETQNS is called in the nonsinglet case and GETQS is called in the singlet 
(or xq + , xF + ) case for calculating the Q 2 evolution. 

5.2 Subroutine GETQNS 

First, the initial distribution is taken from the subroutine QNSO or is read from the 
file #13, and it is stored in the array QNS(I). The distribution is written in the file #12 
if INDIST^l. The initial distribution is stored in the array WQNS(I). The evolved 
distribution at t\ = At is calculated by calling QNSXT and it is stored in QNS(I). The 
major part of the NS evolution calculation is done in the QNSXT subroutine. Next, 
QNS(I) at ti = At is stored in WQNS(I). The distribution at t 2 = 2 At is calculated 
again by calling QNSXT and it is stored in QNS(I). Repeating this step, we obtain the 
final distribution at t — N t At. The results are interpolated either in the variable t or 
in x and they are written in the output file #11. 

5.3 Subroutine GETQS 

The initial distributions, xq s and xg {xqf , xg^, xg A , xqf' A ), are taken from the 
functions, QSO and GO (QIO, QSAO, GAO, QIAO) or from the file #14, #15, #16, or 
#17. The distributions are stored in the arrays, QS(I) and G(I) (QI(I), QSA(I), GA(I), 
QIA(I)), and they are written in the file #12 if INDIST^l. The initial distributions 
are stored in WQS(I) and WG(I) (WQI(I), WDQS(I), WDG(I), WDQI(I), WGHT(I)) 
[see Eqs. (2.20) and (2.21) for xGht and 5(xp ))]. The evolved ones at t\ = At 
are calculated by calling GETQGX, and results are stored in QS(I) and G(I) (QI(I), 
QSA(I), GA(I), QIA(I)). Next, the evolution results are stored in WQS(I) and WG(I) 
(WQI(I), WDSQ(I), WDG(I), WDQI(I), WGHT(I)) and the evolution from h = At 
to t 2 = 2At is calculated. Repeating this step, we obtain the final distributions at 
t = NtAt. The results are interpolated either in the variable t or in x and they are 
written in the output file #11. 

5.4 Subroutines ERR(IERR), INVERT(NN, QQQ) 

The subroutine ERR checks whether the input parameters are valid. If they are not 
valid, this subroutine returns IERR=1 in the main program. The subroutine INVERT 
reverses the array QQQ, namely QQQ(1)->QQQ(NN), QQQ(2)->QQQ(NN-1), 
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QQQ(NN)->QQQ(1). This is used in the Q 2 devolution case (Q 2 < Qg). 

5.5 Functions FUNQ2J(IORDER, ALPHA) and ALPHAN(QQ) 

The function FUNQ2J obtains the Q 2 from a given running coupling constants a s 
within the range from Q 2 - 0.01(GeV 2 ) to Q 2 + 0.01(GeV 2 ) if Q 2 < Q 2 (from Q 2 -0.01 
to Qq + 0.01 if Q 2 < Qq). This function is used in calculating Q 2 from given t. The 
function ALPHAN(QQ) obtains the NLO running coupling constant a s at Q 2 in the 
MS scheme. 

5.6 Subroutine GETQGX(ALPHA) 

The subroutine GETQGX calculates the singlet (and qf , Ff) evolution by calling 
functions, QSXT, GXT, QSKXT, GKXT, GHTXT, DQSXT, and DGXT, which cor- 
responds to the evolution of xq s (xqf), xg, xq s (recomb.), xg(recomb.), xGht, ^(xq s ) 
(5 (xqf)), and 5(xg). (or structure functions) respectively. In the Altarelli-Parisi evo- 
lution in the nucleon, only QSXT and GXT are called. In addition to these, QSKXT, 
GKXT, and GHTXT are called in the Mueller-Qiu case. DQSXT and DGXT are called 
in the nuclear case. 

5.7 Subroutine DELTAOf QGA, QG,DQG) 

This subroutine calculates the difference DQG between the nucleon parton distri- 
bution QG and the nuclear one QGA. 

5.8 Functions QSKXT(I,FF) and GKXT (I) 

These functions calculate recombination effects in the nucleon. QSKXT calculates 
the recombination terms in xq s or xqf (or structure functions), and GKXT calcu- 
lates the recombination in xg. The flavor number FF is used for obtaining either xq s 
(FF=N f ) or xqf (FF=1). 

5.9 Functions QNSXT( I, A LP HA , I ORDER, IT YPE, SIGN), 

QSXT (I, ALPHA, QSK, WWQ,IORS,FF), GXT (I, ALPHA, GK), 
DQSXT (I, A LP HA , WWQ,IORS,FF), DGXT (I, ALPHA) 

These function calculates the Q 2 evolution from tj to tj+At. QNSXT, QSXT, GXT, 
DQSXT, and DGXT calculate the evolution of nonsinglet-quark, singlet-quark, gluon, 
nuclear (actually, nuclear correction of) nonsinglet-quark, nuclear gluon distributions 
(and corresponding structure functions) respectively. The evolution in the brute-force 
method is discussed in detail in section 3. If IORS=l, the singlet distribution is 
calculated, and the xqf (xFf) distribution is calculated if IORS=2. 

5.10 Functions PNS0(I,K,ZK) and PNS1(I,K,ZK,SIGN,WWQ) 

The function PNSO (PNS1) calculates the LO (NLO) nonsinglet splitting function 
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multiplied by a nonsinglet quark distribution or a structure function, which is stored 
in WWQ. If IM0RP=1 (SIGN=-1.0), PNS1 uses the u q - f type splitting function, 
and it uses the u q + f type splitting function if IMORP=2 (SIGN=+1.0). 

5.11 Functions PSQO(I,K,ZK,WWQ,WWG,FF), PSGO(I,K,ZK,WWQ,WWG), 

PSQ1(I,K,ZK, WWQI, WWQ, WWG,FF), and PSG1(I,K,ZK, WWQ, WWG) 

These functions calculate splitting functions multiplied by quark and gluon distri- 
butions. The function PSQO calculates P$(*)WWQ + P$(*)WWG, where WWQ 
is yq s (y), yF s {y), yqt(y), or yFf(y) and WWG is yg(y). The function PSGO calcu- 
late PjP\z)WWQ + P$(2)WWG. PSQ1 and PSG1 calculates the same quantities in 

yy y yy y 

NLO. 

5.12 Functions CQ1(I,K,ZK,ITYPE,WWQ) and CG1(K,ZK,ITYPE,WWG) 

The function CQl calculates the NLO coefficient function multiplied by a quark 
distribution or a structure function (WWQ). CGI calculates multiplied by a gluon 
distribution (WWG). C« and C« are denoted B q n and in Appendix C. 

5.13 Functions EQQ(K,ZK,WWQ), EQG(K,ZK,WWG) and EGQ(K,ZK,WWQ) 

These functions calculate E qq , E qg and E gq multiplied by a structure function or a 
gluon distribution in the one-step evolution. E gg = —E qq is used in calculating the E gg 
part. See Appendix B for the explicit expressions of E qq , E qg and E gq . 

5.14 Function GHTXT(I) 

This function calculates the evolution of the higher dimensional gluon distribution 
xG HT from tj to tj + At in Eq. (2.19d). 

5.15 Subroutines GETFNS(ITYPE, ALPHA, QNS,FNS), 

GETF(IORS,ITYPE, ALPHA, QS, G,FS) 

These subroutines calculate structure functions from quark and gluon distributions 
by taking convolutions with the coefficient functions. These are used in the two-step 
evolution. GETFNS obtains a nonsinglet structure function FNS from a nonsinglet 
quark distribution QNS. GETF obtains a singlet structure function (or xF^~) FS from 
singlet quark (or xqf) and gluon distributions, QS and G. 
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5.16 Function SPENCE(X) 

This function calculates the Spence function by using a series expansion form in 
Eq. (B.8). The upper limit of the summation is taken as max(k)=10000 so that 
accuracy of the obtained function is better than 10 -5 . The Spence function appears 
in the Splitting functions and it could take significant computing time if the function 
SPENCE is repeated called. Therefore, the interpolated one is used in calculating the 
splitting functions. 

5.17 Subroutines D ATARI (N,QFX), DATAR2(N, QFX, GGX), 

DA TAR3(N, QIX, QFX, GGX) 

IF IREAD=2 is chosen, an initial distribution is read from a data file. DATAR1 
reads an initial nonsinglet distribution form the data file #N. DATAR2 reads an initial 
singlet-quark distribution xq s (or structure function xF s ) and a gluon distribution xg. 
DATAR3 reads xqf (F^~), xq s (xF s ), and xg. Then, the initial distribution data are 
interpolated, so that the reasonably large amount of data should be supplied in the 
file. 

5.18 Subroutine SPLINE(N,X,Y,B,C,D) and function SEVAL(N,XX,X,Y,B,C,D) 

The subroutine SPLINE calculates the coefficients, B(I), C(I), and D(I) (I = 1, 2, • • 
•, N) in a cubic Spline interpolation. Y(I) at X(I) (I — 1,2, ■ ■ •, N) are supplied. This 
interpolation program is taken from Ref. |T]J. Using the obtained Spline coefficients, 
the function SEVAL calculates the value of Y at given XX. 

5.19 Functions QQ(X,A,B,C), QQMRS(X,A,B,C,D,E) 

These functions calculate typical parton distributions given by the parameters A, B, 
C, D, and E. QQ calculates Ax B {l-x) c , and QQMRS does Ax B (l+Cy/x~+Dx)(l-x) E . 

5.20 Functions QNSO(X), QSO(X), and G0(X) 

The functions QNSO, QSO, and GO calculate an initial nonsinglet-quark distribution 
xq NS (xF NS ), a singlet-quark distribution xq s (xF s ), and a gluon distribution xg in the 
nucleon. As an example, the MRS(G) distributions [|I(| are provided. The nonsinglet 
one is QNS0=xu„ + xd v =2.70Ax°- 593 (l - 0.76^ + 4.20s) (1 - x) 3 - 96 +0.25 1 3x a335 (l + 
8.63v/i + 0.32x)(l - x) 4 - 41 , the singlet one is QS0=xw t> + xd v + xS=2.704x a593 (l - 
0.76 v ^ + 4.20x)(l-x) 3 - 96 +0.2513x a335 (l + 8.63 v ^+0.32x)(l-x) 4 - 41 +1.74x- 0067 (l- 
3.45v/i+ 10.3x)(l -x) iai , and the gluon distribution is G0=1.51ar a301 (l - 4.14^ + 
10.1x)(l -x) 6 - 06 . 
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5.21 Functions QSAO(X) and GAO(X) 

The functions QSAO and GAO calculate a singlet-quark distribution xq s (xF s ) and 
a gluon distribution xg in a nucleus. As an example, the SK distributions M at 
Q 2 =0.8 GeV 2 are provided. The singlet one is QSA0=1.840x a472 (l - 0.984x) 4 - 06 (l + 
9.33x)+6.423x - 600 (l-x) 813 (l-0.568x), and the gluon distribution is GA0=179.2x L95 (l- 
x) 7 - 32 (l -0.619a;). 

5.22 Functions QIO(X) and QIAO(X) 

The function QIO or QIAO calculates an initial distribution xqf {xF^~) in the nucleon 
or the one in a nucleus. As an example, the MRS(G) and SK xd + distributions are 
provided. The nucleon one is QI0=+0.2513x a335 (l + 8.63v^+0.32a;)(l-x) 4 - 41 +0.4(l- 
0.02)1. 74x-°' 067 (l - 3.45^+ 10.3x) (1 - x) iai +0.043x a3 (l - x) iai (i + 64.9x), and the 
nuclear one is QLA0=1.773x°- 656 (l - 0.956x) 5 - 23 (l + 3.01x)+6.423x a600 (l - x) 8 - 13 (l - 
0.568a;)/3. 
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6. Numerical analysis 



Parton distributions or structure functions either in the nucleon or in a nucleus 
could be evolved in the program BF1. We discuss numerical results on both cases. 

6.1 Accuracy of Q 2 evolution results 

We check the LO and NLO splitting functions, the coefficient functions, and the 
function E by comparing their moments with other calculations in Refs. || 12], |l^ . 



Because the different convention is employed for the gluon distribution from the one in 
Ref. §, and E in Eqs. (2.17) and (2.18) are slightly different. We compare BF1 
evolution results with the evolution of HMRS-B |TJ]]. Choosing the HMRS-B singlet- 
quark and gluon distributions at Qq=4 GeV 2 as the initial distributions, we calculate 
the distributions at Q 2 =200 GeV 2 and compare them with the HMRS-B evolution 
results. We find that both agree well in the sense that the differences are typically 
within 2% in the region (0.01 < x). However, the differences become slightly larger 
at small x (2—4% at 0.0001 < x < 0.001). We also checked various BF1 evolution 
results by comparing with those in the Laguerre polynomial method ||. It indicates 
that both results agree within about 1% accuracy. Furthermore, we checked that 
results in the one-step evolution [xF(x, Qq)—>xF(x, Q 2 )} agree well with those in the 
two-step evolution [xq(x, Qq) — >xq(x, Q 2 )—>xF(x, Q 2 )]- From these comparisons and 
repeated checks on the program, the program BF1 is believed to be a reliable evolution 
program. We discuss the details of numerical accuracy in the following. 

First, Q 2 evolution of a nonsinglet distribution is calculated in the program BF1. 
There are essentially two parameters which determine numerical accuracy of the Q 2 
evolution. These are numbers of points in variables x and t (N x , N t ). We show 
dependency of our results on these numbers. For running the program BF1, the input 
parameters and an input distribution in section 4 should be supplied. As an initial 
nonsinglet distribution, we employ the HMRS-B xu v + xud at <3o=4 GeV 2 . It is evolved 
up to Q 2 =200 GeV 2 with A=0.19 GeV and JV)=4. 

In Fig. la, N x is fixed at ^=1000, then jV t =10, 20, or 500 is taken for finding 
the dependency on N t . For example, the input parameters in the case of A^=500 are 
IOUT=l, INPUT=1, IREAD=1, MODIFY=l, INDIST=2, IORDER=2, ITYPE=4, 
ISTEP=1, IMORP=l, Q02=4.0, Q2=200.0, DLAM=0.255, NF=4, XX=0.0, NX=1000, 
NT=500, NSTEP=100, NXMIN=-4, NA=1, RFM=0.0, and DKHT=0.0. There are 
solid, dashed, and dotted curves at Q 2= 200 GeV 2 ; however, these cannot be distin- 
guished in Fig. 1. It indicates that numerical results are almost independent of the 
number Nt in the nonsinglet case if Nt is more than a certain number. This result is, in 
fact, anticipated because the scaling hypothesis works approximately in parton distri- 
butions. In comparison with the A^=500 results, numerical differences are within 1% 
at 10" 4 < x < 0.5 and 4% at x=0.8 in the ^=10 case, within 0.6% at 10~ 4 < x < 0.5 
and 2% at x=0.8 in JV,=20, and within 0.2% at 10" 4 < x < 0.5 and 0.8% at £=0.8 in 
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N t =50. 

Next, N x is varied with the fixed N t =50 in Fig. lb. Four types of curves are 
shown at Q 2 =200 GeV 2 . The dotted curve is the evolution result for N x =100, the 
dashed for jV x =300, the dot-dashed for A^=1000, and the solid for JV t =4000. As it is 
shown in the figure, we should take rather larger number of N x for getting converging 
results. The accuracy is worse if N x <100. In comparison with the N x =A000 results, 
numerical differences are within 9% at 10~ 4 < x < 0.5 and 24% at x=0.8 in the 7^=100 
case, within 3% at 10~ 4 < x < 0.5 and 8% at x=0.8 in JV X =300, and within 0.7% at 
10" 4 < x < 0.5 and 2% at x=0.8 in JV X =1000. 

From the above results, we recommend to use N t >50 and N x >1000 for getting 
accuracy better than 1% at 0.0001 < x < 0.5 and better than 2% at 0.5 < x < 0.8 
in the nonsinglet case. Running CPU time with N t =50 and A^=1000 is about five 
minutes on the AlphaServer 2100 4/200. 

Singlet distribution results are shown in Fig. 2a, where N t is varied with fixed 
iV^lOOO. The initial distributions are the HMRS-B singlet-quark and gluon distribu- 
tions at Qo=4 GeV 2 . The dotted, dashed, dot-dashed, and solid curves are the results 
for iV t =20, 50, 200, and 1000 respectively. As it is obvious from the figure, we need 
to take large N t for getting convergent results at small x. In comparison with the 
iV(=1000 results, numerical differences are within 5% at 10~ 4 < x < 0.9 in the N t =20 
case, within 2% at 10" 4 < x < 0.9 in At=50, and within 0.4% at 10" 4 < x < 0.9 
in N t =200. In Fig. 2b, N x is varied with fixed N t =200 in the singlet distribution. 
The dotted, dashed, dot-dashed, and solid curves are obtained by taking N x =100, 300, 
1000, and 4000 respectively. In comparison with the ^=4000 results, numerical differ- 
ences are within 8% at 10~ 4 < x < 0.5 and 24% at x=0.8 in the case ^=100, within 
3% at 10" 4 < x < 0.5 and 8% at x=0.8 in A^=300, and within 0.7% at 10~ 4 < x < 0.5 
and 2% at x=0.8 in A^=1000. 

Evolution results of a gluon distribution are shown in Figs. 3a and 3b, and they 
show similar tendency to those in Figs. 2a and 2b. In. Fig 3a, Nt is varied with fixed 
^=1000. In comparison with the ^=1000 results in Fig. 3a, numerical differences are 
within 5% at 10~ 4 < x < 0.6 in the N t =20 case, within 2% at 10" 4 < x < 0.6 in N t =50, 
and within 0.4% at 10~ 4 < x < 0.6 in JV t =200. In Fig. 3b, N x is varied with fixed 
N t =200. In comparison with the 7X^=4000 results in Fig. 3b, numerical differences 
are within 25% at 10~ 4 < x < 0.6 in the case A^=100, within 8% at 10~ 4 < x < 0.6 
in A^=300, and within 1% at 10" 4 < x < 0.2 and within 2% at 0.2 < x < 0.6 in 
^=1000. The accuracy is slightly worse in the gluon distribution but it is mainly in 
the medium x region, which is not important in the gluon distribution. 

From these results, we recommend to use N t >200 and N x >1000 in the singlet 
case for getting accuracy better than 1% at 0.0001 < x < 0.5 and 2% at 0.5 < x < 0.8 
in the singlet-quark distribution, and accuracy is better than 1% at 0.0001 < x < 0.2 
and 2% at 0.2 < x < 0.6 in the gluon distribution. Typical running time with N t =200 
and ^=1000 on the AlphaServer is about one hour in the NLO singlet case. 

We also check numerical accuracy of our evolution results in the case of % + % 
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distributions, in the case with the recombinations, and in the case of nuclear parton 
distributions. However, the results have similar dependency on N t and N x , so that 
they are not discussed in this paper. A user can always check one's evolution results 
by changing N t and N x if one is suspicious about the accuracy. 

The one-step evolution, e.g. xq s (x, Ql)^>xq s (x, Q 2 ), is used for getting the above re- 
sults. We check the accuracy in the two-step evolution case, e.g. xq s (x, Ql)^xq s (x, Q 2 ) 
-^xF s (x, Q 2 ), by calculating the F 2 structure function of the proton with the MRS(G) 
input distributions. We find that the accracy is better than 1% in the region (10~ 4 < 
x < 0.8) by taking A^=200 and 7^=1000. Similar results are obtained in the recombi- 
nation case, so that we do not discuss their numerical details. 

From the above discussions, we had better use N t >200 and N x >1000 for getting 
accuracy better than 2%. Typical running time could be an hour on the AlphaServer. 
If the program is run on a SUN workstation, it could take a significant amount of 
running time. Therefore, a user may use a reasonably powerful machine in running the 
program BF1. 



6.2 Comparisons with experimental data 

We first compare nonsinglet evolution results with CDHSW -F 3 -structure-function 



data |15 in neutrino reactions. We choose the initial distribution as the MRS(G) 
xu v +xd v distribution at Q 2 =A GeV 2 flHJ, XM 1 ,+xc/ t) =2.704x a593 (l-0.76 v ^+4.20x)(l- 
x) 3 - 96 +0.2513x a335 (l+8.63v/i+0.32x)(l-x) 4 - 41 , The two-step part of the evolution pro- 
gram is used for calculating xF^(x, Q 2 ) from the initial xu v (x, Ql)+xd v (x, Qq). Chosen 
input parameters are the following in the NLO case: IOUT=2, INPUT=1, IREAD=1, 
MODIFY=l, INDIST=1, IORDER=2, ITYPE=3, ISTEP=2, IMORP=l, Q02=4.0, 
Q2=400.0, DLAM=0.255, NF=4, XX=0.045, NX=1000, NT=50, NSTEP=100, NXMIN 
=-4, NA=1, RFM=0.0, and DKHT=0.0. The LO evolution is also calculated by 
choosing IORDER=l. The results are shown in Fig. 4, in which the solid (dashed) 
curves indicates the NLO (LO) calculations at x=0.045, 0.225, and 0.55. The LO 
results are different from the RK-MK-SK results in Ref. || because the initial distri- 
bution is different. The LO and NLO xF% are already different at the starting point 
(Q 2 =4 GeV 2 ) due to the NLO contributions from the coefficient function. The figure 
shows that the BF1 evolution agrees with the experimental Q 2 dependence. As it is 
obvious from the figure, the NLO contributions become significant at small Q 2 (~1 
GeV 2 ). 

Next, evolution of the proton F 2 structure function is studied. As it is shown 
in the previous subsection, the BF1 program can be run at small x, as small as 
x = 10~ 4 , with good accuracy, so that evolution results could be compared with re- 
cent HERA F 2 data. The initial distributions are the MRS(G) distributions at Qo=4 
GeV 2 . First, we calculate F 2 ^ d +(x, Q 2 ) and F 2t s{x,Q 2 ) by setting the initial ones as 
QI0=+0.2513a; - 335 (l + 8.63 v ^+0.32a;)(l-x) 4 - 41 +0.4(l -0.02)1. 74x- a067 (l-3.45v^+ 
10.3x)(l-x) iai +0.043x a3 (l-x) iai (l+64.9x), xg s =2.704x a593 (l-0.76v^+4.20x)(l- 
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x) 3 - 96 +0.25 1 3x a335 (l + 8.63 v ^ + 0.32x)(l-x) 4 - 41 +1.74x-°- 067 (l-3.45 v ^+10.3x)(l- 
x) 101 , and xg=1.51a; _0 ' 301 (l — 4.14y / a7 + 10.1x)(l — x) 6m with the input parameters: 
I0UT=6, INPUT=1, IREAD=1, M0DIFY=1, INDIST=1, I0RDER=2, ITYPE=2, 
ISTEP=2, IM0RP=2, Q02=4.0, Q2=200.0, DLAM=0.255, NF=4, XX=0.0013, NX 
= 1000, NT=200, NSTEP=100, NXMIN=-4, NA=1, RFM=0.0, and DKHT=0.0. The 
scale parameter is A=0.255 GeV and the number of flavor is four in the evolution. Next, 
F 2 , s +(x, Q 2 ) is calculated with xs+ =0.2(1 - 0.02)1. 74a;-°- 067 (l - 3.45^+ 10.3z)(l - 
x) 10A . These evolution results are added to give Ff (x, Q 2 ) = — (l/3)[F 2>d +(x, Q 2 )+F 2>s + 
(x,Q 2 )] + (4:/9)F2 t s(x,Q 2 ). The LO Altarelli-Parisi and Mueller-Qiu evolution results 
are also obtained by changing IORDER and MODIFY. These results are shown in Fig. 



5 with various Ff data at x ^0.0013, 0.013, 0.13, and 0.45 from SLAC M, 111 (x=0.14, 



0.45), BCDMS m (x=0.U, 0.45), EMC pj, (x=0.125, 0.45), NMC [|(|(x=0.0125, 
0.14), Fermilab-E665 [g| (x=0.0012, 0.012), ZEUS J22| (x=0.0016, 0.014, 0.11), and 
HI |23j (x=0.0013, 0.013, 0.13). The solid, dashed, dot-dashed curves are obtained in 
the LO Altarelli-Parisi, NLO Altarelli-Parisi, and NLO Mueller-Qiu evolution equa- 
tions respectively. The NLO Q 2 dependence agrees with the experimental data as 
shown in the figure, except for the Fermilab-E665 at small x and at small Q 2 where 
perturbative QCD would not work. At this stage, it seems that Ff data are not ac- 
curate enough to find the recombination effects on the evolution. We do not step into 
the details of recombination analysis in comparison with the HERA data. 

Finally, Q 2 evolution of a nuclear structure function F 2 A is studied. Q 2 depen- 
dence of the ratio R^ a = F^/F® had been measured by NMC f2Hfl . We com- 
pare our evolution with the NMC data in Figs. 6a, 6b, 6c, and 6d at x=0.0085, 
0.035, 0.125, and 0.45. We assume that F 2 D is equal to the nucleon one for sim- 
plicity. The initial distributions at <5o=0.8 GeV 2 in the nucleon are taken from 
Ref. [§: xq NS = xu v + x4=i-967x a501 (l - x ) 389 (l + 9.27a;), xg s =1.967x a501 (l - 
a;) 3 - 89 (l + 9.27a;)+2.058a; - 294 (l - a;) 1L2 (l + 7.95a;), and xc/=16.11:r - 839 (l - a;) 5 - 29 (l - 
0.597a;). Those in the calcium nucleus are xg^=1.840x 0472 (l - 0.984x) 406 (l + 9.33x), 
xgf a =1.840x°- 472 (l - 0.984x) 4 - 06 (l + 9.33a;)+6.423a; ' 600 (l - a;) 8 - 13 (l - 0.568x), and 
xg Ca =l79 .2a; 1 ' 95 '(1 — a;) 7 ' 32 (l — 0.619a;). Both nonsinglet and singlet structure func- 
tions are calculated. Then, we obtain F 2 A =(l/18)F 2 A 7V5 +(2/9)F 2 j4 s for isoscalar nu- 
cleus because SU(3)fi aV or symmetry is assumed in antiquark distributions in Ref. 0]. 
In Figs. 6a, 6b, 6c, and 6d, the solid, dashed, dot-dashed curves are obtained in the 
LO Altarelli-Parisi, NLO Altarelli-Parisi, and NLO Mueller-Qiu evolution equations 
respectively with A=0.2 GeV and Nf=3. As shown by the figures, NLO and recom- 
bination contributions to the ratio are conspicuous at small x (=0.0085, 0.035). They 
are very small at medium x (=0.45) in Fig. 6d; however, it does not mean that their 
contributions to the structure functions themselves are very small. If we evolve F 2 from 
Qq=0.8 GeV 2 , the recombination effects are larger than the NLO ones. It is interesting 
to find such large recombination contributions in Fig. 6a. However, the recombination 
cannot be tested at this stage because we do not have the data in the wide Q 2 region 
at small x. 
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7. Conclusions 



We investigated numerical solution of Q 2 evolution equations of spin-independent 
parton distributions (and structure functions) in the nucleon and in nuclei by using the 
brute-force method. Two types of evolution equations are studied. One is the usual 
Altarelli-Parisi equations and the other is the modified ones due to parton recombina- 
tions (Mueller-Qiu). 

We divide the variables x and t into small steps and simply solve the evolution 
equations. Numerical accuracy depends essentially on the numbers N x and N t . We 
obtain excellent numerical solution in the wide x range (0.0001 < x < 0.8) by taking 
large enough N x and N t . We recommend to use N x > 1000 and N t > 50 in the 
nonsinglet case and N x > 1000 and N t > 200 in the singlet for getting numerical 
accuracy better than 2% in the x range, 0.0001 < x < 0.8. 

We provide the FORTRAN program BF1, which can be run by supplying the input 
parameters and the input distribution(s). Parton distributions xq NS , xq s , xqi + xqi 
(corresponding structure functions), xg, and those in nuclei could be evolved in the 
program. The program can also handle the Q 2 devolution from Qq to smaller Q 2 . 
This is a very useful program in studying structure functions in the nucleon and in 
nuclei. Typical running CPU time is five minutes on the AlphaServer 2100 4/200 in 
the nonsinglet and one hour in the singlet. Therefore, a reasonably powerful machine 
should be used for running the BF1 program in the singlet case. 
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Appendix A. Running coupling constants 

Running coupling constant in the leading order (LO) is 



* L S °{Q 2 " 



Air 



/? ln(Q 2 /A 2 ) 

and the one in the next-to-leading order (NLO) is |26[ 



(A.l) 



«f LO (Q 2 ) 



47T 



/3 ln(Q 2 /A 2 ) 



fa lnln(Q 2 /A 2 ^ 



(A.2) 



In Eq. (A.2), the renormalization scheme is MS and A is the QCD scale parameter in 
this scheme. (3q and (3i are given by 



(3 = ^C G -~T R N f , (5 x = ^C 2 G -^C G N f -2C F N f , 



with the color constants 

C G = N C 



N 2 c -l 
2N r 



N c is the number of color (N c =3) and Nf is the number of flavor. 



(A3) 



(AA) 



Appendix B. Splitting functions 

Splitting functions in the leading order are 



4 0) (*) 


= c F 


1 + x 2 
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= 2Tr 


' x 2 + (1 




^ l + (l-x) 
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= 2C G 
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5(1 -x) 



+ - — - + x(l - x) + 
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where the + function is defined by 



11 1 N f T R 

12 3 C G 



(B.l) 



6(1 -x) 



o (1 — X) 



/o 1 — X 

It should be noted that the above integration is defined in the region < x < 1 



(B.2) 
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The NLO splitting functions for qf(x) and g(x) are given by Herrod and Wada 



P 



<? + (*) = Sa P^ NS (x) + 2C F T R F qq (x) 



P ( l } (x) 



- C G T R F l (x) - C F T R F 2 (x) , 



P {1 Ux) 



91i 



-(%F%(x) - C F C G F 2 q (x) - C F T R N f F* q (x) 
C G T R N f F l gg {x) + C F T R N f F 2 gg {x) + C G F*(x) 



<// // -■ • 7.7' • ■• 7.7- (B.3) 

, T 1 

77'-- ' - ■■■ (>-'./■■ 77'- ■ - 'J^gg^ 
The NLO nonsinglet splitting function for a "q+g type" distribution [q NS = ^ aj((fc ± <&)] 

i 

is given by f| 

"3 1 
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+ -C f C a ^P g (x)tPa(x) + 5(1-x) 



2 tt 2 + C(3)-8S(oo) 
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where P F (x), P G (x), P/v F (x), and Pa(^) are given by Curci, Furmanski, and Petronzio 
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The ( function is defined by ((k) = — - and the numerical value (£(3)=1. 2020569...) 

n=l n 

is taken from Ref. p7| . S'(oo) is given by the ( function as S(oo) = — §C(3). The 
integral in Pa{x) is expressed in terms of the Spence function S(x): 
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It should be noted that other convention is sometimes used, namely —S(x) may be 
called the Spence function |[27|| . It is useful to use a series expansion form for numerical 
calculations |27[| : 

S(x) = - £ • (B-S) 

k=i K 



The functions F qq , F qg , F gq , and F 99 in Eq. (B.3) are defined by 
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4 „„ o 20 
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where I x is the Spence function [I x = S(x)] and J x is defined by 

J x = -— 7r 2 -lnxln(l+x) + 1 S(l + x) . (5.10) 

In the one-step evolution equation in Eq. (2.14), the matrix E is introduced. Each 
matrix element is given by 

E qq (x) = -E( lgq C 9 ) , 

E qg {x) = E{ lqq C°)-E{ lqg Ci)-E{ lgg C°) , (Bn) 

E gq (x) = E{ lgq C«) , 

E gg (x) = -E qq (x) . 



where 



E( lM C)/2C F T R = -| + ^-|-^+(l + 5x-lx 2 )lnx 

+ (1 + x) (in 2 x - 21nx In (1 - x) + 21^ 
+ (-^-l+x + ^x 2 )ln(l-x) , 
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A and i? in the above equation are 
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In the above equation, we corrected the factor (17/3 + 170x/3 — 40x 2 /3), which was 
(17 + 170x/3 - 40x 2 /3) in tdie original paper |§. 
The splitting function P qg = xP qg is given by 



2 (-2x + 15x 2 - 30x^ + 18x 4 



(5.14) 
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Appendix C. Coefficient functions 



The coefficient functions C% 9 (x, a s ) in Eqs. (2.10) and (2.11) are written by the 
functions BS,' 9 (x): 



C*(x, a.) = 5(1 -x) + -±B*Jx) , C<l(x, a.) = ^B«(x) 



(CI) 



The quark part is given by 



C h 



Bf(x) = -±[F q {x)-4x] , 

B q 2 {x) = C F F q (x) , 
Bl{x) = C F [F q (x) -2-2x] 



where the function F g (x) is given by 

F q (x) = 



2(l-x)+ 
The gluon part is: 



3 1 + X 2 1 , , 1 + X 1 , 

f-(9+5x)-2y— - lnx+2(l+x 2 ) 



ln(l — x) 



1 — x 



Bf(x) = 2T R [F g (x) - 4x(l - x)} 
B 9 2 {x) = AT R F g (x) , 



where 



F g (x) = (1 -2x + 2x 2 )ln 



1 — x 



x 



(C.2) 



-5(1-x)(9+|tt 2 ) 



2 

3' 
(C.3) 



+ 8x(l - re) - 1 . 



(C.4) 
(C.5) 
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Figure Captions 



Fig. 1 (a) N t dependence of nonsinglet evolution results is shown. Next-to-leading- 
order Altarelli-Parisi evolution results are calculated by running the program 
BF1. JV X =1000 is fixed and N t is varied (A t =10, 20, and 500) in order to check 
numerical accuracy due to the choice of N t . The initial distribution is the HMRS- 
B xu v + xd v distribution at Qq=4 GeV 2 . There are dotted, dashed, and solid 
curves at Q 2 =200 GeV 2 for N t =10, 20, and 500 respectively; however, they cannot 
be distinguished in the figure because differences among them are fairly small. 
See text for the details of the input parameters, (b) N x dependence is shown by 
taking fixed N t =50. The dotted, dashed, dot-dashed, and solid curves are the 
results for JV X =100, 300, 1000, and 4000 respectively. 

Fig. 2 (a) N t dependence of singlet-quark evolution results is shown. Next-to-leading- 
order Altarelli-Parisi evolution results are calculated. A^ x =1000 is fixed and N t 
is varied. The initial distributions are the HMRS-B xq s and xg distributions at 
<5o=4 GeV 2 . The dotted, dashed, dot-dashed, and solid curves at Q 2 =200 GeV 2 
for iV t =20, 50, 200, and 500 respectively, (b) N x dependence is shown by taking 
fixed iV t =200. The dotted, dashed, dot-dashed, and solid curves are the results 
for A^=100, 300, 1000, and 4000 respectively. 

Fig. 3 (a) N t dependence of gluon evolution results is shown. N x is fixed at A^=1000. 
The input parameters and the input distributions are the same with those in Fig. 
2. The dotted, dashed, dot-dashed, and solid curves at Q 2 =200 GeV 2 are for 
N t =20, 50, 200, and 500. (b) N x dependence is shown by taking fixed A t =200. 
The dotted, dashed, dot-dashed, and solid curves are the results for 7X^=100, 
300, 1000, and 4000 respectively. 

Fig. 4 Nonsinglet Q 2 evolution results are compared with CDHSW data. The initial 
distribution is the MRS(G) xu v + xd v at Qo=4 GeV 2 . The solid and dashed 
curves are the NLO and LO results at x=0.045, 0.225, and 0.55. 

Fig. 5 Q 2 evolution results of the proton's F 2 structure function at x=0.0013, 0.013, 
0.13, and 0.45 are compared with various data in SLAC, BCDMS, EMC, NMC, 
Fermilab-E665, ZEUS, and HI experiments. The solid, dotted, and dashed curves 
are obtained by using the NLO Altarelli-Parisi, LO Altarelli-Parisi, and NLO 
Mueller-Qiu equations respectively. F$ at £=0.0013, 0.013, and 0.13 are shown 
by 10 • F%(x = 0.0013), 10 2/3 • F%(x = 0.013), and 10 1/3 • F%(x = 0.13). 

Fig. 6 Q 2 dependence of the structure function ratio F2 a (x,Q 2 )/F®(x,Q 2 ) is calcu- 
lated and results are compared with NMC data. The solid, dotted, and dashed 
curves are obtained by using the NLO Altarelli-Parisi, LO Altarelli-Parisi, and 
NLO Mueller-Qiu equations (A=0.2 GeV, N f =3) respectively at (a) x=0.0085, 
(b) 0.035, (c) 0.125, and (d) 0.45. The input distributions in the nucleon and in 
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the calcium nucleus are taken from Ref. [|J. For example, AP (NLO) means that 
the next-to-leading-order Altarelli-Parisi equations are used in evolving F 2 in the 
nucleon and the one in the calcium, then the ratios are taken. 
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TEST RUN OUTPUT 



IOUT= 3 INPUT= 1 IRE AD = 1 MODIFY = 1 INDIST= 1 

IORDER= 2 ITYPE= 4 ISTEP= 1 IMORP= 2 

Q02= 4.0000 Q2= 200.000 DLAM= 0.2550 NF= 4 

XX= 0.0000000 NX=1000 NT= 200 NSTEP= 50 NXMIN= -4 

NA= 1 RFM= 0.00 DKHT= 0.000 



0.000100 
0.000120 
0.000145 
0.000174 
0.000209 
0.000251 
0.000302 
0.000363 
0.000437 
0.000525 
0.000631 
0.000759 
0.000912 
0.001096 
0.001318 
0.001585 
0.001905 
0.002291 
0.002754 
0.003311 
0.003981 
0.004786 
0.005754 
0.006918 
0.008318 
0.010000 
0.012023 
0.014454 
0.017378 
0.020893 
0.025119 
0.030200 
0.036308 



18.335139 
17.128091 
15.996575 
14.936139 
13.942574 
13.011907 
12.140385 
11.324464 
10.560801 
9.846243 
9.177817 
8.552722 
7.968325 
7.422150 
6.911875 
6.435325 
5.990469 
5.575416 
5.188411 
4.827834 
4.492194 
4.180129 
3.890402 
3.621892 
3.373589 
3.144582 
2.934037 
2.741170 
2.565211 
2.405346 
2.260649 
2.129987 
2.011924 



66.652696 
61.663113 
57.000809 
52.646335 
48.581338 
44.788502 
41.251490 
37.954894 
34.884181 
32.025646 
29.366368 
26.894167 
24.597561 
22.465735 
20.488499 
18.656257 
16.959975 
15.391153 
13.941795 
12.604385 
11.371860 
10.237586 
9.195334 
8.239259 
7.363868 
6.564002 
5.834797 
5.171655 
4.570205 
4.026253 
3.535733 
3.094648 
2.699013 
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0.043652 
0.052481 
0.063096 
0.075858 
0.091201 
0.109648 
0.131826 
0.158489 
0.190546 
0.229087 
0.275423 
0.331131 
0.398107 
0.478630 
0.575440 
0.691831 
0.831764 
1.000000 



1.904595 
1.805604 
1.711949 
1.620025 
1.525750 
1.424866 
1.313447 
1.188551 
1.048903 
0.895387 
0.731157 
0.561477 
0.393883 
0.239332 
0.113441 
0.033130 
0.002843 
0.000000 



2.344797 
2.027887 
1.744079 
1.489125 
1.258859 
1.049428 
0.857641 
0.681417 
0.520250 
0.375536 
0.250490 
0.149345 
0.075655 
0.029963 
0.007968 
0.001021 
0.000023 
0.000000 




Figure 1 (a) 




Figure 1 (b) 
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Figure 3 (a) 




Figure 3 (b) 
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